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Abstract 

We consider a probabilistic cellular automaton to analyze the stochas- 
tic dynamics of a predator-prey system. The local rules are Markovian 
and are based in the Lotka-Volterra model. The individuals of each 
species reside on the sites of a lattice and interact with an unsymmet- 
rical neighborhood. We look for the effect of the space anisotropy in 
the characterization of the oscillations of the species population den- 
sities. Our study of the probabilistic cellular automaton is based on 
simple and pair mean-field approximations and explicitly takes into 
account spatial anisotropy. 
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1 Introduction 



In the last years a particular great effort has been done in order to understand 
the role of space given by a spatial structure and local interactions in the 
characterization of the dynamics of competing biological species systems [U 

I21IHISIISIIS1IZIIHIIS1 HOI HH OS IXSl OS OZZl IXB|- in this context it 

has been studied irreversible stochastic lattice models [HI (20J [21] with the 
purpose of mimic predator-prey systems with Markovian local rules based 
in the Lotka-Volterra model J22[ [23]. One of the problems that has been 
object of study is the connection between the time oscillations of population 
densities and spatial pattern distribution of the individuals of each species. 

Here we study the coexistence of a two-species system by considering a 
probabilistic cellular automaton (PCA) which is a modified version of the 
automaton devised in [TUJ [HI [T71 [18] . The model, to be called anisotropic 
predator-prey PCA possess local rules that are similar to the ones of the 
cellular automaton proposed in [21] and was introduced by us in order to 
explore the effect of spatial anisotropy in the temporal oscillations. 

We report dynamic mean-field approximations which take into account 
the spatial dependence of densities and pair correlations of sites. In the 
next section we present the model. In Sec. 2 we show the equations for the 
time evolution of species densities and show the spatial dependence of these 
equations. In Sec. 3 and 4 we perform dynamic simple and pair mean-field 
approximations. Last section summarize the model, method and results. 

2 Model 

The physical space occupied by the species is represented by a regular lattice 
of iV sites in which each site can be in one of three states. At each site of 
the lattice we attach a stochastic variable rj i that takes the values 0, 1 and 
2 according whether the site i is empty, or occupied by a prey or occupied 
by a predator, respectively. The state of the system can be represented by 
the vector rj = (rj^ rj 2 , ■ ■ ■ , rj N ). The time evolution equation for Pi(rj), the 
probability of state rj at time £, is given by 

P t+ \r ] )=Y J W(r } \r l ')P\r l ') (1) 

a' 

where the sum is over the 3^ configurations of the system and W(rj\r)') is the 
transition probability from a state rf to state ry, given that at the previous 
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Figure 1: Transitions of the predator-prey model. The three states are: prey 
(1), predator (2) and empty (0). The allowed transitions obey the cyclic 
order shown. 



time step the system was in state if. Since we are considering probabilistic 
cellular automata, all the sites are updated simultaneously. In this case we 
have 

N 

W( V \v') = l[w l ( Vl \ V ') (2) 

8=1 

where Wii^^rf) is the conditional transition probability per site. 



2.1 The anisotropic predator-prey PC A 

The stochastic rules, embodied in the transition rate u^t^I?/), are set up 
in order that the allowed transitions between states are only the ones that 
obey the cyclic order shown in Fig. [H Prey can only born in empty sites; 
prey can give place to a predator in a process where the prey dies and the 
predator is instantaneously born; finally predator can die leaving an empty 
site. The empty sites are places where prey can proliferate and can be seen 
as the resource for prey surveillance. The death of predators complete this 
cycle, reintegrating to the system the resources for prey. 

The anisotropic predator-prey PCA has three parameters: a, the proba- 
bility of birth of prey, b, the probability of birth of predator and death of prey, 
and c, the probability of predator death. Two of the process are catalytic: 
the occupancy of a site by prey or by a predator is conditioned, respectively, 
to the existence of prey or predator in the neighborhood of the site. The 
third reaction, where predator dies, is spontaneous, that is, it occurs, with 
probability c, independently of the neighbors of the site. We assume that 
a + b + c — 1 with < a, b, c < 1. 

The transition probabilities of the anisotropic predator-prey PCA are 
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described in what follows: 

(a) If a site % is empty, 77^ = 0, and if at north or east there is at least one 
prey, then it can be occupied by a prey, r q i = 1, in the next time step, with 
a probability proportional to the parameter a and to the number of prey n a 
at north and east of site i. 

(b) If a site is occupied by a prey, 7^ = 1, then the site has a probability 
of being occupied by a new predator, % = 2, in the next time step if there 
are prey at north or east. In this process the prey dies instantaneously. The 
transition probability is proportional to the parameter b and the number of 
predators at north and east of site %. 

(c) If site i is occupied by a predator, 7^ = 2, it dies spontaneously with 
probability c. 

The anisotropic cellular automaton is a variation of the automaton intro- 
duced in [TUl [HI [T71 HH] . Here each site of a regular square lattice interacts 
with its first neighbors only at two preferential directions. This anisotropic 
neighborhood consists of the northern and eastern neighbors of each site as 
shown in Fig. [2J The set of transition probabilities per site is given by 

Wl (l\ V ) = ^mtyfc.O) + (1 - ln 2 )8(ViA), (3) 

w i {2\ v ) = ~n 2 6( Vi ,l) + (l-c)S( Vi ,2) (4) 

and 

Wl (0\v) = (1 - \a ni ) 6(7],, 0) + c%, 2) (5) 

where 

n i = ^2 S (Vk^) and n 2 = ^2$(Vk> 2 ) ( 6 ) 

fc fc 

and the sum is over the neighbor sites localized at east and north of site i and 
correspond to the number of neighbors of site i occupied by prey individuals 
and predators individuals, respectively. 

We have considered this probabilistic cellular automaton with the purpose 
of verifying the effect of anisotropy in the properties of the time oscillations 
of the predator-prey system. The rules considered are in some sense inspired 
in the north-east-center (NEC) cellular automaton [22], which also consider 
an unsymmetrical neighborhood of northern and eastern sites. 

The present stochastic dynamics predicts the existence of states, called 
absorbing states, in which the system becomes trapped. Once the system 
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Figure 2: A site (C) of the square lattice and its two nearest neighbor sites 
to the east (E) and to the north (N). The layers, n, n+l, n + 2, ... are 
perpendicular to the southwest-northeast direction. 



has entered such a state it cannot escape from it anymore remaining there 
forever. There are two absorbing states. One of them is the empty lattice. 
The other absorbing state is the one in which the lattice is full of prey. This 
situation occurs if there are few predators and they become extinct. The 
remaining prey will then reproduce without predation filling up the whole 
lattice. The existence of absorbing stationary states is an evidence of the 
irreversible character of the model or, in other words, of the lack of detailed 
balance [2T]. However, the most interesting states, the ones that we are 
concerned with in the present study, are the active states characterized by 
the coexistence of prey and predators. 

2.2 Evolution equation for state functions 

The densities of prey, predator and empty sites are defined as 

Pt(l) = (5(T,i,l)h (7) 
Pf(2) = (5( Vi ,2)) e , (8) 

Ff(Q)= (8( Vi ,Q))t, (9) 

respectively. The lower index i is used to denote the site and the upper index 
£ stands for the time. The pair correlation of two neighbor sites % and j one 
being occupied by prey and the other by predator is defined as 

4(1 2 ) = (6( Vi ,l)S( Vj , 2))* (10) 
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In this definition the sites i and j are two neighboring horizontal sites, the 
site j being at the east of i. For two neighboring sites placed vertically we 
use the notation 



where k denotes the neighbor of i to the north. 

The evolution equation for the densities can be obtained by using the 
rules of the automaton and the evolution equation for the probability P l (rj). 
Using equations (j2J), (J3J and we obtain the time evolution equation 
for the density of prey at site i, 



where we have used prime and unprimed to denote quantities at time I + 1 

and i, respectively. Again, P^-( ) denotes a pair correlation of a site i which 

is empty and its neighbor at the north j which is occupied by a prey and 
Pik( 1 ) denotes the pair correlation of the site i which is is empty and 
its neighbor to the east k which is occupied by a prey. Note that site i is the 
site to be updated. The time evolution equation for the density of predators 
at site i is given by 



The evolution equations for correlations of two neighbor sites are given 
by equations which involves clusters of two, three and four neighbors to the 
north and east and they are too cumbersome. We will write them in the pair 
approximation in the next subsection. 

We will assume in the present analysis that the densities and the pair cor- 
relations are homogeneous so that the correlations of sites placed horizontally 
or vertically are the same. In other words the systems exhibits a specular 
symmetry along the southwest-northeast line (see Fig. [2]). Therefore, 




(11) 




(12) 




(13) 
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(14) 



and 




(15) 



6 



3 Simple mean-field approximation 



We will first analyze the equations (|T2|) and (fT3l) for the densities of prey 
and predators by means of a simple mean field approximation [2TI [26| 127] . 
In this approximation we write the probability of a given cluster of sites as 
the product of the probabilities of each site. In our analysis we maintain the 
space dependence of each probability. This is necessary since the rules that 
define the model are not isotropic. We will use the notations 

x n = P n (l), y n = P n {2) z n = P n {0). (16) 

where the index n stands for a site at the n layer. A layer is defined as the 
sites belonging to a line perpendicular to the southwest-northeast axis, as 
shown in Fig. [2j Using this approximation and considering equations ffT2l) 
and f|T3|) we obtain 

x' n = az n x n +i + x n — bx n y n +i, (17) 

and 

y' n = bx n y n+ i + (1 - c)y n . (18) 

where z n = l — x n — y n . Due to this relation there is no need for the equation 
for the density of empty sites. 

The analysis of the above set of equations show that the stable solutions 
are of two types: a prey absorbing state, where the all lattice is full of prey; 
and an active solution with prey and predators densities constant in time 
and space. So, we have just homogeneous and nonoscillating solutions when 
we treat the anisotropic PCA by means of simple mean-field approximation. 
At this level of mean-field approximation the space anisotropy does not play 
any role in determining the kind of coexistence of species. 

4 Pair- mean field approximation 

Now we consider the evolution equations for one-site correlations and two 
site (pair) correlations. The evolution equations for the pair correlations 
depend on clusters of three and four sites. Following the approach of the 
pair approximation [25], 12S1 [27J [28] we write the clusters of three and four 
sites as products of pair correlations and one-site correlations. However, 
we must be careful in doing this approximation because we are maintaining 
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the spatial anisotropy dependence of the one and two site correlations. For 
example, the correlation of a site at layer n in state 0, a site at layer n + 1 in 
state and a site at layer n + 2 in state 1, is written in the pair approximation 

as 

p /nmA _ ^^+1(00)^+1^+2(01) nQ x 
^n,n+l,n+2(UUlJ - — . [LU) 

There are more complex clusters that appear in the evolution equation for 
pair correlations. For example, the following cluster of four sites is approxi- 
mated by 

t-> Z \ r n,n+l 

(02)P nin+1 (00)P n+1 , n+2 (01) 

^n.n+l.n+2 ^ Q 1 ^ ~ P(0)P n+1 (0) ' K ' 

Observe that P„ n+ i(10) and P„ n+ i(01) are two independent variables. 
Now we use the notation 



Pn,n+l(01) 




Pn,n+l(12) 


= v n , 


Pn,n+l(02) 


= w n , 


(21) 


Pn,„+l(10) 


= fni 


Pn,n+l(21) 


= 9n, 


Pi,n+l(20) 


Jim 


(22) 


Pn,n+l(H) 




Pi,n+l(22) 


Sri) 


Pn,n+l(00) 


= q n - 


(23) 



Using the pair approximation we derive the following equations for pair 
correlations, 

/ / QnU n fnQn u n+l\ , /-, a \ I , u n v n+l 

u n = a — a— J + I 1 - «) ( «n - b— 

V z n z n z n+l J z \ x n+l 

^h n u n+ i 

2 \ Z n Z n X n +i J z n+l 

+ c(g n - b 9 -^) . (24) 

V X n+1 J 

, ab f U n V n+1 f n U n Vn+l\ . 1 n JnW„ 
V n = IT I 1 J + 9 a (l - C ) 
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•Efi+i 2 x. 



n+l 



+ (l-c)((l--K-6-2ti), (25) 
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Figure 3: Spatial oscillations of prey and predator along the southwest- 
northeast direction. 



w n = b[ (i-o)-!— - - 27ir~: 

u 2 n v n+1 ( au n w n+1 

-a h (1 - c) \w n - - 



2 Z n 



+ 6c^±i + c(l- C ) Sn , 

u n q n u n q n u n+1 \ ac u n w r , 
a- 1 



fnUn b ( V n f n V n f n U n +i \ 



^n+i 2 \ rr 



a- 



n+1 / 



+ c I (1 - j)w„ - j — 



, _ abv n f n u n+1 l ft / t> n r» _ b v n r n v n+1 



n 



+o(l - c)- 



(26) 



(27) 



9 



:i-c)(g n -b 9 -^), (28) 

b ( v nfn v nfn u n+l 



H + (1 - c) ft n - a 



+ c(l - c)s n . (29) 
The equations for the densities are 

x^ = au n + x n -bv nJ , (30) 



and 



2/n = + (1 - c)y n , (31) 



Therefore, we get a set of eight equations that we have iterated to find 
the solutions. We considered periodic boundary conditions and an initial 
condition where half of the lattice with a set of densities and the other half 
with another set of densities. 

For a = b and great values of c the system attains a stationary state which 
is the absorbing prey state. Decreasing c there is a transition to an active 
state, coexistence of species, and we found very interesting solutions which 
are characterized by travelling waves of the densities of each species. In Fig. 
[3]we show, for a particular set of the parameters a = b and c = 0.05, the prey 
and predators oscillations as a function of space and for a fixed instant of 
time. We can see that the behavior of the space oscillations are very complex 
and very different from a sinusoidal wave. 

We remark that our assumption concerning the spatial dependence of 
densities and pair correlations underlying our mean-field approach could not 
easily be conceived a priori. It was set up by considering the type of lo- 
cal dynamics with unsymmetrical rules and from the results of Monte Carlo 
simulations for the present PCA [10, 29J. To further clarify this point we 
show in Fig. H] a snapshot generated by simulation of the present PCA on 
a square lattice where each site can be occupied by prey, predators or can 
be empty. Periodic boundary conditions were used and the system evolved 
in time according to the rules defined by Eqs. Q, (TjJ and (jSJ). A syn- 
chronous update was used. We see that the distribution of individuals of 
each species is not homogeneous but exhibits a pattern composed by layers 
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Figure 4: Snapshot of a configuration obtained from Monte Carlo simulation 
of the PCA in a square lattice with periodic boundary conditions. Predator 
and prey individual are represented by black and grey points, respectively. 
The fronts move from northeast to southwest. 

of prey and predators. These layers are displayed in a manner perpendicular 
to the southwest-northeast direction, the direction along which the oscil- 
lations occur as can be seen in Fig. HI We have found that the spatial 
oscillations do not occur independently of the time oscillations, since the 
spatial layers of prey and predators are not static but are like fronts moving 
along the southwest-northeast direction. The spatial pattern oscillations are 
intimately associated to local time oscillations of the species pm EH] . These 
features were incorporated in the present mean-field analysis of the model. 

5 Conclusions 

We have considered a predator-prey probabilistic cellular automaton with 
anisotropic local stochastic rules. We studied this model by means of dy- 
namic mean-field approximation at two orders: simple mean-field approx- 
imation and pair mean-field approximation. Due to the anisotropy these 
approximations only can be performed if we consider the spatial dependence 
of probabilities. The simple mean-field approximation for this automaton 
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just provides the prey absorbing state and active spatial homogeneous so- 
lutions which are also constant in time. The pair mean-field approximation 
gives much more rich results and show that the active states characterized 
by time oscillations of species densities are inhomogeneous in space. There- 
fore, we have spatiotemporal patterns of coexistence. This is in accordance 
with previous Monte Carlo simulations [TU] where we have found local time 
oscillations connected to inhomogeneous spatial distributions of species indi- 
viduals. 
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